if (!require("tidyverse")) install.packages("tidyverse")
if (!require("rms")) install.packages("rms")
if (!require("rmda")) install.packages("rmda")
if (!require("ggDCA")) devtools::install_github("yikeshu0611/ggDCA")
if (!require("ggplot2")) install.packages("ggplot2")
if (!require("ranger")) install.packages("ranger")

rm(list = ls())
df_select <- readRDS("~/analysis/lyz_ml/rds/step_03_mlr3_select.RDS") %>% as.data.frame()

colnames(df_select) <- make.names(names(df_select))
df_select$group <- as.numeric(df_select$group) - 1

row_sample <- sample(nrow(df_select), round(nrow(df_select)*0.75))

df_train <- df_select[row_sample,]
df_test <- df_select[-row_sample,]

# logistic 与 随机森林 DCA数据获取
dca_fun <- function(train,test, group = "group", type = "logistic") {
  train <- train %>% dplyr::select(group, everything())
  test  <- test %>% dplyr::select(group, everything())
  colnames(train)[1] <- "group"

  
  if (type == "logistic") { # 1 Logistics 回归
    model <- glm(group ~ ., binomial(link='logit'), train )
    test$logistic <- predict(model, test[, -1], type = "response") %>% as.numeric()
    df_model <- ggscidca::dca(test, # 指定数据集,必须是data.frame类型
      outcome = "group", # 指定结果变量
      predictors = "logistic", # 指定预测变量
      probability = T,
      graph = F
    )
  }
  else if (type != "logistic") { # 2 随机森林
    model <- ranger(group ~ ., train)
    test$randomForest <- predict(model, test[, -1], type = "response")$predictions
    df_model <- ggscidca::dca(test, # 指定数据集,必须是data.frame类型
      outcome = "group", # 指定结果变量
      predictors = "randomForest", # 指定预测变量
      probability = T,
      graph = F
    )
  }
  
  df_model$net.benefit %>%
    pivot_longer(
      cols = -threshold,
      names_to = "model",
      values_to = "NB"
    ) %>%
    return()
}

# 合并数据
df_model <- rbind(dca_fun(train = df_train,test = df_test,group = "group",type = "logistic"), 
      dca_fun(train = df_train,test = df_test,group = "group",type = "randomForest") %>% subset(model == "randomForest"))

# df_model <- rbind(dca_fun(train = df_select,test = df_select,group = "group",type = "logistic"), 
#                   dca_fun(train = df_select,test = df_select,group = "group",type = "randomForest") %>% subset(model == "randomForest"))

# 重新排序
df_model$model <- factor(df_model$model, levels = c( "randomForest","logistic", "all", "none"))
df_model %>% ggplot(aes(threshold, NB, color = model)) +
  geom_line(linewidth = 1.2) +
  scale_y_continuous(limits = c(-0.02, 1), name = "Net Benefit") +
  scale_x_continuous(limits = c(0, 1), name = "Threshold Probility") +
  theme_classic(base_size = 16) +
  # theme(
  #   legend.position = c(0.2, 0.3),
  #   legend.background = element_blank()
  # ) +
  ggsci::scale_color_jama(name = "Models")

ggsave(filename = "result/step_07_dca_test.pdf", width = 32, height = 20, dpi = 300, units = "cm")

# ggsave(filename = "result/step_07_dca_all.pdf", width = 32, height = 20, dpi = 300, units = "cm")

